Wigner crystallization in a single and bilayer graphene 
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We study the possibility of Wigner crystallization in both single- and and bi-layer graphene using 
a real space tight binding model. In addition to verifying our earlier prediction for single layer 
graphene, we predict that the bilayer graphene can undergo a Wigner crystal transition at low 
enough carrier density thus offering a possibility of tunable and bipolar Wigner crystal. We find 
that aside from the Coulomb interaction of the carriers of two layers that stabilizes the charge 
ordered state, so does the inter layer coupling. 

PACS numbers: 



The study of low dimensional system covers a great 
deal of area in the research field of physics. In the 
past couple of years, the experimental success of fab- 
ricating single [1, 2] and bi-layer graphene [3] has en- 
riched the interest of studying the properties of low di- 
mensional strongly interacting electronic systems both 
theoretically and experimentally. The low dimensional 
system of graphene is more interesting because many of 
the properties of these systems are quite different from 
that of the conventional two dimensional system; such as 
low temperature [3, 4] and room temperature [5] quan- 
tum Hall effect, and suppression of weak localization of 
the carriers [6, 7]. From technological point of view, 
theses systems open up possibilities of atomic-scale elec- 
tronic devices which can go beyond the limitations of the 
silicon-based electronics. In addition to the small size of 
these materials, the very high mobility of the carriers of 
the single and bilayer graphene looks to be promising for 
devices which depend up on the transport properties of 
the carriers. 

The importance of studying the single- and bi-layer 
graphene not only comes from the fact that they have 
different physical properties compared to the other con- 
ventional system such as two dimensional electron gas 
(2DEG) prepared in the heterostructure of GaAs and In- 
GaAs but also from the fact that their properties are very 
different between themselves. One of the best example 
is the experimental verification of the different quantum 
Hall effect in the single [4] and bilayer graphene [3] . Both 
single and bilayer system show weak localization effect 
but with different minimum zero-bias conductivity. The 
theoretical prediction of the Klein paradox in the single 
layer graphene is another example [8]. 

Structurally, graphene has a hexagonal lattice struc- 
ture having two inequivalent lattice sites, A and B, per 
unit cell. Bilayer graphene is obtained from two iden- 
tical single layer graphene with a z-axis stacking in a 
specific way known as Bernal stacking [9]. The intra 
layer nearest-neighbor hopping energy is about ten times 
stronger than the inter layer hopping. So, to the lowest 



order approximation, each layer of the bilayer graphene 
is almost similar to the single layer graphene, except that 
the two layers are coupled via small inter layer hopping. 
The effective difference can be seen in the low energy 
dispersion relation of the charge carriers in the conduc- 
tion band which for the single and bilayer graphene can 
be written as, £k = ifiv/k, and Ek ~ ±t_L ± y/t^ + £fc 
,[9] respectively, where, we set h = 1. These dispersion 
relations are obtained using tight binding model where 
yp = ^ = 5.SeVA is the Fermi velocity, t = 3.2eV 
is the intra layer nearest-neighbor hopping energy, and 
tj_ ~ OAt is the inter layer hopping energy. For bilayer 
graphene for small k, the dispersion relation for two of 
the branches can be approximated as, Ek ~ ^2^- ^'^ 
due to the inter layer hopping the massless particles of 
the single layer graphene become massive where the mass 
of the carriers, m, is related to inter-layer hopping en- 
ergy by, m — t±/vp. An important point to note here is 
that the dispersion relation of the single layer graphene 
is quite different than that of the 2DEG system whereas 
the dispersion relation of the bilayer graphene has similar 
form, i.e., parabolic dispersion, as that of 2DEG system. 
Flat dispersion in bilayer effectively lowers the kinetic en- 
ergy for particles with small momenta. This effect will 
become important for our discussion of Wigner Crystal- 
lization (WC). 

Here we want to study the charge ordering effects in 
single versus bilayer graphene. Difference in electronic 
spectra could manifest itself as a difference in the possi- 
ble charge ordered states. Indeed we find this to be the 
case. We find that WC can occur in a bilayer system at 
small doping limit in zero magnetic field. On the other 
hand the single layer graphene is stable against charge 
ordering [10]. WC is a physical phenomenon envisioned 
by Wigner [11], in which the crystallized phase appears 
when electrons localize and form a crystal to minimize 
the potential energy, while paying the concomitant ki- 
netic energy cost which arises from localization. This 
phase arises in 2DEG system which has parabolic dis- 
persion relation as we decrease the carrier density, n, of 
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FIG. 1: Lattice structure of graphene single layer having 60 
lattice sites that we used to solve the hamiltonian. The index 
number of some of the lattice sites are shown. 



the system because the ratio of the potential to kinetic 
energy varies as which can be made as large as it is 
needed to get the crystallized phase by making the sys- 
tem more and more dilute. 

This discussion brings obvious questions about the 
WC in graphene: a) Can a single layer graphene have 
a Wigner crystal phase even if it has linear dispersion 
relation? b) Is it obvious that the bilayer graphene will 
have Wigner crystal phase since it has parabolic disper- 
sion realation? c) If bilayer can have crystallized phase 
what is the structure of the crystallized phase? We will 
answer these questions by calculating the occupation (n) , 
of carriers on lattice sites of single and bilayer graphene 
using a real space tight binding model. If {n) is equal for 
all sites, we will define the state as homogeneous phase, 
and if it not equal for all sites, we define the state as 
some sorts of inhomogeneous phase, one of which is the 
Wigner crystal phase. We will determine the structure 
of the inhomogeneous phase with varying (n). 

We now use a real space tight-binding model a) to ver- 
ify that single layer graphene remains in the liquid phase, 
and b) to predict that bilayer graphene can undergo a 
crystallization transition. We study the electronic distri- 
bution on the lattice of single layer and bilayer graphene 
having 60 atomic sites in each layer. The real space dis- 
tribution of the system of 60 lattice sites that we have 
chosen and some of the corresponding indices is shown 
in Fig. 1. In addition to this layer, the bilayer graphene 
has another layer below or above this layer with a shift of 
the whole layer by one lattice spacing in the lateral direc- 
tion. For the second layer we index the lattice sites with 
indices 61 to 120. We use the real space tight binding 
hamiltonian. 

The tight-binding Hamiltonian to be solved is as fol- 
lows: 

k—l,i,j ky^l.i.j 

^ (1) 

ki^lj i ki 



where i,j denote the lattice sites in both single and bi- 
layer, and fc, Z = 1, 2 represent index for the first or sec- 
ond layer of bilayer graphene. For single layer, k ~ I = 1. 
Since the first term of Hamiltonian has summation over 
«,j only for k — I, tkijj represents the hopping energy 
of the electrons only within a layer. We do calculation 
including only the nearest neighbor hopping with energy 
equal to tki.ij = t = 3.2eV for all lattice sites. For the 
second term of the Hamiltonian the summation is only 
for k I which implies that for single layer graphene this 
term is zero and for bilayer graphene t±Aki,Aij represents 
the inter layer hopping energy. In this term A denotes 
one of the non-equivalent lattice sites. Since the bilayer 
graphene can be obtained by the special stacking called 
the Bernal stacking of two single layers grpahene, elec- 
trons can hop between two layers only through one of the 
sublattices which we have chosen to be A. For the sum- 
mation of this term, since one of the two layers is shifted 
laterally by one lattice spacing, if i is an even numbered 
lattice index, j will be an odd numbered lattice index or 
vice-versa. Vq is the on site coulomb repulsion, Vki,ij is 
the interacting potential between two electrons at posi- 
tion rki and rij , is the chemical potential, and Uki , nij 
are the occupation number of the electrons on lattice sites 
i,j of the layer k,l respectively. Except otherwise indi- 
cated, we use Vki n = ^ exp(-|rfc,-ri, |/2a) .^]^gj,g 2a is the 

screening length, so that we could regularize the on-site 
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potential Vq — Vki=ij = to be finite. There are indi- 
cations that Coulomb terms are decaying faster than 1 /r 
that would also be consistent with our approximation of 
taking into account only a few lattice sites [12] [13]. 

We use mean-field theory to simplify the Coulomb in- 
teraction term. The Hamiltonian then can be expressed 
as, 

HmF = - ^ tki,ljcliCij - ^ t±Aki,Al]C'Af^iCAlj 
k—l,i,j k^l,i,j 

+ ^{Wki + Uki ~ li)nki - C, 

(2) 

where Wu = EyVfe,j,K), Uk^ = C = 

lJ2k:^ijVkt,ij{nkt}{nij} + EH^o("HT)W»i) is a con- 
stant. 

Using periodic boundary condition we write down a 
Hamiltonian matrix corresponding to, 

hki,lj — -~tki,lj — t±Aki,Alj + {Wki + Uki — fJ^)Ski,lj, (3) 

for single and bilayer graphene. For single layer graphene, 
we calculate Wki including carriers up to next-next near- 
est neighbor. For 60 lattice sites in a single layer, the 
Hamiltonian becomes a 60 x 60 symmetric matrix. For 
the bilayer graphene also, we calculate the contribution 
to Wki coming from both the layers fc, I including the 
next-next nearest neighbor. The Hamiltonian for the bi- 
layer graphene becomes a 120 x 120 symmetric matrix. 
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FIG. 2: (Color online) Charge distribution in ordered state of 
bilayer graphene is shown. Larger circles correspond to larger 
charge density on lattice sites. The sites connected by solid 
lines belong to the top layer and by dotted lines belong to 
bottom layer. The left figure, which we define as A-phase, 
has charge ordering where larger charge densities are sitting 
on sub-lattice A. The right figure, which we define as B-phase, 
has charge ordering where large charge densities are sitting on 
sub-lattice B. 

To find the distribution of electrons on the lattice we 
start with a random distribution of {uki). Using random 
initial input of {riki), we construct h]^i^ij and find the 
eigenvalues £'„ and eigenvectors (jfj.^. We calculate new 
{nki) using following relation, 

n 

En 

where, f{En) = l/(e''B^ -I- 1) is the Fermi distribution 
function. Using the standard procedure of the mixing of 
the old and new (rifci), we repeat the calculation until we 
reach the convergence of {uki)- We check that the results 
we obtain do not depend on initial distribution of riki- 

a) Single layer case. We calculate occupation num- 
ber, (nj), of electrons on each lattice site of the single 
layer graphene. We choose different values of the chemi- 
cal potential and Coulomb potential to vary the average 
occupation number below half filling and evaluate {rii) 
for all i. We always find rii = nj for all sites. This de- 
fines a homogeneous charge distribution also known as 
the liquid phase. Thus we verify our earlier prediction 
that the single layer graphene always remains in the liq- 
uid phase. Here we remark that this homogeneous phase 
assumes that graphene sheet has no ripples [14] . 

b) Bilayer case. Similar to the case of single layer 
graphene, we calculate (ui) for each lattice site. In bi- 
layer case, in addition to the Coulomb energy and chemi- 
cal potential, we can tune the inter layer hopping energy, 
t± . First we imagine turning off the Coulomb energy and 
study the charge distribution as a function of tj_ . We find 
that (rii) — (rij) for all j — i + 2, and (rii) ^ {nj) for all 
j = i + 1. This implies that all A-sites have equal occu- 
pation. Similarly, all B-sites have equal occupation. But, 
site A and B do not have equal occupation. This defines 
a commensurate triangular inhomogeneous phase. Since 
this inhomogeneity is driven only by the inter-layer hop- 
ping energy, we define this phase as kinetic energy driven 



(KED) inhomogeneous phase. Depending up on the aver- 
age occupation of the charge carriers on the lattice, also 
known as filling factor, we find two distinct inhomoge- 
neous phases, KED-A and B phase. A schematic repre- 
sentation of the charge distribution in these two phases is 
shown in Fig. 2. We find KED-A phase when the filling 
factor is less than half filling and the KED-B phase when 
it is more than half filling. We find a cross over between 
these phases at half filling. So exactly at half filling the 
system will be in liquid phase. 

When we turn off t± and turn on Coulomb energy, we 
again find the inhomogeneity in the charge distribution. 
The charge distribution resembles to that of the B-phase 
as shown in Fig. 2. Since this phase is driven by the 
inter-layer Coulomb interaction, we define this phase as a 
Coulomb interaction driven Wigner crystal (WC) phase. 
The pattern of the inhomogeneous charge distribution 
does not change below and above the half filling case; of 
course the system regains the liquid phase far above the 
half filling. 

If we include both the inter layer Coulomb interaction 
and the inter layer hopping, we find competing phases 
below the half filling and cooperating phases above it. 
Above the half filling KED-B and WC-B phases cooper- 
ate with each other. Below the half filling, KED-A and 
WC-B phases compete with each other. Since this com- 
peting phase is more interesting, we are going to explore 
it more for a system which has 40% filling. 

To define the degree of the inhomogeneity of both the 
A and B phases we define a charge order parameter, A = 
'^^riB+nA^ . In principle this order parameter is position 
dependent, but in our calculation it is transactionally 
invariant. For B-phase it gives a positive number and for 
the A-phase it gives a negative number. We will be using 
the absolute value. The result of the calculation of this 
order parameter as functions of t± and Vq is shown in 
Fig. 3. 

For a finite Vq, an increase in t± decreases the order 
parameter of WC-B phase. A further increase in t± de- 
stroys the WC-B phase turning it in to a liquid phase and 
finally drives the system to a KED-A phase. Similarly, 
for a finite t± for zero Vq the system is in the KED-A 
phase. Increase in Vq destroys the KED-A phase turn- 
ing it to a liquid phase and finally drives the system in 
to WC-B phase. The black colored region on the figure 
separates the two competing phases. For realistic param- 
eters we expect Vq ^ t± and WC-B state always wins. In 
practice disorder and fluctuations will produce mixture 
of these two states in real systems. 

How one can detect WC in a bilayer? One exper- 
iment that would be sensitive to the charge ordering 
at atomic scale is scanning tunneling spectroscopy. It 
would allow one to measure spatial variations of local 
density of states. To better characterize WC states we 
calculate the local density of states (LDOS) Nk{ri,uj) = 
\'f'ki^i)\'^^i^ ~ -^n), in the WC phase for bilayer and 
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FIG. 3: (Color online) The order parameter of competing 
KED-A and WC-B phases as functions of t± and Vb at (n) — 
0.4. Increasing Vo increases the order parameter of WC-B 
phase, whereas increasing t± increases the order parameter of 
KED-A phase. The black-colored region represents the phase 
boundary between theses two competing phases. 



in liquid phase for a single layer graphene, see Fig. 4. 

DC transport measurements is a less direct alternative 
way to detect WC phase. Depinning effects and other 
transport anomalies could be another way to detect sug- 
gested WC phase. 

In conclusion, we show that the single layer graphene 
always remains in a liquid phase. We also show that the 
bilayer graphene is bound to undergo Wigner crystalliza- 
tion at reasonable value of the Coulomb interaction. We 
find that in principle there are two crystallized configura- 
tions possible in a bilayer. The in-phase KED-A state is 
stabilized by interlayer hopping. The out-of-phase WC- 
B state is stabilized by the Coulomb interactions. LDOS 
maps, that can be measured with STM, would reveal 
charge alternations between sites and would allow direct 
imaging of the WC state in bilayer graphene. 
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